Contour integral method for thermal and quantal fluctuations 
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The partition function by means of the static path approximation (SPA) plus the random-phase 
approximation (RPA) treatment can be written as a contour integral form without solving the RPA 
equations for a separable interaction. This method is an efficient way to evaluate numerically the 
partition function for realistic calculations. As an illustration, we adopt the pairing model at finite 
temperature. 
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Mean-field approximations, such as the finite tem- 
erature Hartree-Fock Bogolyubov (FTHFB) approach 
J, 0, 13 , are standard methods for describing nuclei at 
finite temperature, and the collective motions around 
the static HFB solution are obtained by solving the fi- 
nite temperature random phase approximation (RPA). 
At zero temperature, the HF plus RPA method gives a 
good approximation to ground-state energy of the exact 
shell model calculation Q. However, the FTHFB-I-RPA 
method is not applicable to systems at high temperature 
because of large thermal fluctuations. The inclusion of 
fluctuations and correlations beyond thermal mean field 
can be done by use of the path integral representation 
of the partition function. The static path approxima- 
tion (SPA) H S is useful treatment to evaluate ap- 
proximately the partition function in finite systems with 
separable interactions, and provides the exact result at 
high temperature. However, as temperature decreases 
the SPA becomes inaccurate because the quantal fluctu- 
ations around the mean fields can no longer be neglected. 
Small amplitude fluctuations around such a configuration 
give corrections to the partition function similar to the 
standard RPA around the mean field. The SPA+RPA 
method 0, 0, can take the effects into account. 
It should, however, be pointed out that in the limit of 
zero temperature the formalism breaks down since it in- 
cludes contributions of small-amplitude quantal fluctua- 
tions around unstable mean-fields. 

In realistic calculations for deformed nuclei or heavy 
nuclei, the evaluation of the RPA partition function is 
difficult due to time-consuming computation, because a 
complete set of RPA frequencies and amplitudes must be 
known with high accuracy for each temperature and the 
number of the RPA eigenvalues is very large. Recently, 
a contour integral representation for the standard RPA 
correlation energy |l2i | and more efficient way {13i] based 
on the linear response theory have been proposed. These 
methods are valuable and useful for evaluating the stan- 
dard RPA correlation energy. In this letter, we present a 
contour integral method for the RPA partition function 
at finite temperature. 

We consider the total Hamiltonian 



where Hq and V are the single-particle Hamiltonian 
and a two-body residual interaction, respectively. Here, 
the two-body interaction is assumed to be of separable 
multipole-multipole form interaction: 
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where Qa is a one-body Hermitian operator 

Qo. = (3) 

kl 

Using the Hubbard-Stratonovich transformation : the 
functional integral representation for the partition func- 
tion is written as an auxiliary field path-integral 

Z = Tr{e-^") 

= j P[(T„(r)]exp ^-i^xa^ cr2(T)dT^ 



cTr[fexp(- / drH^r)], 



H'{t) = Hq -^XaO-aiT)Qa, 
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where (3 — 1/T is the inverse temperature, and T de- 
notes time ordering operator. Let us now use the Fourier 
transformation of the auxiliary fields CTq, defined as 



O-air) = + Vane' 
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with the Matsubara frequencies uj„ — 2'im/j3 and the 
static variables CTq,. The SPA is obtained by consider- 
ing only the static paths in the evolution operator. The 
Hamiltonian H'{t) can be divided into the static and 
fiuctuation parts: 



H'{t) = Hs+5H{t), 



(7) 



where Hg and SH{t) are the static and fluctuation part 
as follows: 



H = Ho + V, 



(1) 



5H{t) = - Xa VanQge' 
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Substituting Eqs. (6-9) into Eq. (4) and performing the 
Gauss integral over the coefficients rjan, we obtain 
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Crpa = Y[det{l-xR{iOn))-\ 
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where £fe arc the cigcnenergies of H', and -R(w„) is the 
response function matrix given as 

I? \ - \" ^kllkl^kljfk - fl) i^^s 
iiaa'(i^n) - 2^ 2 , , ,2 ' (1^) 
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with £ki = £k — £; and the Fermi occupation probabil- 



ity /fe = (1 + e"^'') . Solving the dispersion equation 
det(l — xR{iuJn)) = or diagonalizing the finite temper- 
ature RPA equations, Crpa can be written exactly in 
terms of the RPA frequencies f^j/, 
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RPA = 



-J-sinh%i 
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(13) 



where the prime in Yi'ia restricts the product to pairs 
{k,l) that satisfy k < I and Ski 0. Note that for de- 
formed nuclei or heavy nuclei there are many RPA fre- 
quencies rii^ in the numerator and s^i in the denominator. 
For arbitrary a, the lowest Q^^'s may become zero and 
imaginary (or complex) at low temperature. However, 
Crpa{o') remains finite and positive both for Q,^ = 
and for imaginary Qu if f3\fl\ < 2tt. This condition will 
be violated at very low temperature. 
Let us define the thermal RPA energy 



Erpa = — 



dlnCppA 
d(3 



= ^ ( ^■'Coth^ - ^ Efejcoth^ j , 

\ I' k<l / 

(14) 

where we neglected the temperature dependence of the 
RPA frequencies Civ In the limit (3 oo, the thermal 
RPA energy Erpa reduces correctly the RPA correlation 

energy 

ERPAiT = 0) = i(^Q,-^£,,J. (15) 

\ 1^ k<l J 

We now introduce the analytical function of the com- 
plex vaiable z 



f{z) = det(l - xRiiz)) 
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FIG. 1: Contours used in evaluating the integrals (19) and 
(20). Crosses denote the positions of all the RPA roots fl^ 
and the poles £ki 



As mentioned above, the FTRPA frequencies fl^ are ob- 
tained as the zeros of the determinant, det(l — X-R(iw)) = 
0, i.e., f{u) = 0. Here, we notice that u> = ±Sli, are the 
roots of /{lu) = and the function f{uj) has obviously 
poles at a; = ±eki ■ Thus we obtain the following relation 



= #log/(^) 



dz 
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with a derivative f{zy = df{z)/dz. Applying Cauchy's 
theorem to a meromorphic function f{z), we obtain the 
following relation 
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dzg{2 



f{z) 



Y.9{eki)-Y.9{^")^ (18) 



k<l 



where g{z) is an arbitrary complex function which is ana- 
lytical in the enclosed region. The contour of the integra- 
tion goes aromid all positive roots and all poles of f{z) 
as shown in Figs. 1(a) and (b), where we consider the 
case of real RPA frequencies. The thermal RPA energy 
Erpa, Eq. (14), is written as the closed contour integral 
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= ^Jdzizcoth^)^, 



(19) 
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(20) 



where the closed path for Eq. (19) is taken so that the 
function zcothf3z/2 is analytical in the closed region(see 
Fig. 1(a)). A modified contour integral for Eq. (20) is 
illustrated in Fig. 1(b), where the contribution from the 
semicircle vanishes in the limit of infinite radius. These 
are the extended form of the contour integral of the RPA 
correlation energy at zero temperature 12, 13]. Integrat- 
ing Erpa = —dlnCfiPA/dp over (3 in Eq. (14) provides 
the RPA partition function Crpa 
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In comparison with the usual RPA method of Eq. (14), 
the advantage of the contour integral form in Eqs. (19- 
20) is that one can choose the integration path so that 
the integrand is smooth enough and the mesh of numer- 
ical integration along such path can be much larger than 
the actual energy intervals of the FTRPA solutions Oj^. 
The contour integral calculation reduces drastically the 
computation time without loss of precision. In practical 
calculations of the above method, we evaluated the con- 
tour integral (20) illustrated in Fig. 1(b) neglecting the 
function {zcol\vj3z /2)' . The calculated results in Figs. 2- 
4 agree well with those of the SPA plus RPA calculation 
with Eq. (14). However, in general the neglect of the 
factor (zcoth/3z/2)' is not justified in all cases. In that 
case, this factor should be suitably taken into account in 
Eqs. (19) and (20). 

As an illustration, we consider a monopole pairing 
Hamiltonian 



H = 



k 



(22) 



with k the time reversed states. Here, Sk is a single- 
particle energy and P is the pairing operator P — 
J2k(^kCk- By means of the SPA+RPA |ll Inlbased 
on the Hubbard-Stratonovich transformation [lj|) tbe 
canonical partition function is given by 



Tr[PAre-^/^]sPA+RPA 



(23) 
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^fcSinh[Afc/r] 
2AfeSinh[Ljfe/2r] ' 



(24) 
(25) 



where we introduced the number parity projection Pjq = 
(1 -|- (Te*'^^)/2 (cr means the even or odd number 
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FIG. 2: (Color online) Thermal energy and thermal corre- 
lation energy as a function /3. The dashed line denotes the 
SPA, the soUd line the SPA-I-CIR, and the broken line the 
SPA-hRPA. 



parity) [T^ instead of the exact number projection. Here 
use the notation ojk for the conventional thermal RPA 
energies and Aj, — ^ e'^ + , eJc ffc - - G/2, 
and 7fe — £fe — /i — Afc. The SPA partition function 
is obtained by neglecting the RPA partition function 
Crpa- Then the thermal energy can be calculated from 
E = —d\nZ'ff/df3. In this calculation, we use the single- 
particle energies Ek given by an axially deformed Woods- 
Saxon potential with spin-orbit interaction |15| . The 
Woods-Saxon parameters are chosen so as to approxi- 
mately reproduce the experimental single-particle ener- 
gies extracted from the energy levels of the odd nucleus 
^^■^Sn (^^^Sn core plus one neutron). The deformation 
parameter estimated from the experimental B{E2) value 
is (5 = 0.23 in the even-even nucleus ^^-^W. The 50 doubly 
degenerate single-particle levels are taken outside ^^^Sn 
core, and we fix the pairing force strength at G = 20/^ 
MeV. As mentioned above, the SPA-f-RPA breaks down 
at low temperature. However, it has recently been shown 
that in the monopole pairing case the SPA+RPA with 
the number parity projection reproduces well exact re- 
sults for low temperature jinl |. Thus, the number parity 
projection is essential to describe the thermal properties 
for low temperature. 

For the partition function (23), we performed the 
three cases of numerical calculations; the SPA, the SPA 
plus contour integral representation (SPA+CIR) with 
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FIG. 3: (Color online) Heat capacity as a function of (i. The 
open squares denote the SPA, the solid circles the SPA+CIR, 
and the open circles the SPA+RPA. 



FIG. 4: (Color online) Level density as a function of thermal 
energy E. The same symbols as those in Fig. 3 are used. 



Eq. (19-20), and the SPA plus RPA (SPA+RPA) with 
Eq. (14) using the RPA frequencies. The thermal ener- 
gies are calculated from E = —dh\Z'^/d(3. The thermal 
energy obtained in the SPA+CIR (SPA+RPA) deviates 
from that of the SPA. The energy difference between the 
SPA+CIR (SPA+RPA) and the SPA is the correlation 
energy Ecorr ^ E — Espa- Figure 2 shows the thermal 
energy E and the RPA correlation energy Ecorr as a func- 
tion of p. We can see that the SPA+CIR results agree 
well with the SPA+RPA one. Thus, the contour integral 
method is useful for studying the thermal and quantal 
fluctuations. For comparison, we solved the Richardson 
equations and calculated exactly the ground-state 

energy Egg. The obtained energy Egg = -74.48 MeV is 
quite close to the SPA+CIR value, E = -75.24 MeV at 
/3 = 10 MeV"\ supporting the SPA+CIR method. 

We can see that the heat capacity exhibits the char- 
acteristic S-shape behavior as shown in Fig. 3. This 
S-shape heat capacity was recently observed [T^ [T9l | in 
^^^Dy, ^^^Er, and ^^^Yb, and was interpreted as an sig- 
nature of the breaking of nucleon Cooper pairs. The 
SPA+RPA and SPA+CIR results deviate from the SPA 
one in higher temperature region above T =0.5 MeV, 
and the quantal fluctuations are important for this re- 
gion. In our previous paper [23, |^ , we clarified that the 
S-shape behavior of heat capacity is attributed to the 
reduction of the pairing energy which can be calculated 
from G(Ptp) = GTd\nZ%/dG. 



When problems are associated with the variable iV, 
they can be solved by considering a canonical ensemble, 
in which the particle number is strictly fixed. Since we 
used the number parity projection instead of the exact 
number projection, the level density for a system with N 
particles is given by an inverse Laplace transformation of 
the partition function in Eq. (23) 
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p{E,N) = — / 



00 + ioa 



d(3eP^Z%{(3). (26) 



In the saddle point approximation for the f3 integral, the 
level density is given by 



p{E,N) 



al3E 
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(27) 



Figure^ldepicts the natural logarithm of the level density 
(|27|l as a function of the thermal energy. We can see that 
the SPA+CIR result agrees with the SPA+RPA one. 

In conclusion, we developed the contour integral 
method which drastically reduces the computation time 
without loss of precision when evaluating the thermal 
and quantal fluctuations. This method is the extension 
of contour integral representation at zero temperature 
to the thermal system, and is numerically efficient. In 
this paper, we considered the simple pairing model as an 
illustration. The contour integral method is also applica- 
ble to more realistic interactions such as the pairing plus 
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quadrupole-quadrupole interaction. This study is now in progress. 
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